## Note: please run Study2_1_structural_topic_models.R prior to this file and save STM_result_X.Rdata.

#### setting environment ####
require(quanteda)
require(stm)
require(stringi)
require(runjags)
require(coda)
require(psych)

newspaper.names <- c("Asahi", "Chugoku", "Chunichi", "Hokkaido", "Kahoku", 
                     "Mainichi", "Nikkei", "Nishinippon", "Sankei", "Yomiuri")

#### data preparation ####
# document-feature matrix of newspaper editorials published from October 1, 2017 to September 30, 2018
load("Study2_dfm.Rdata")
# estimation results of the factor analysis models
load("FA_result.Rdata")

#### estimates by the 65-topic model ####
# point estimates (posterior means)
Study2.ideal.points.mean <- rowMeans(t(rbind(FA.result$mcmc[[1]][, 1:10], 
                                             FA.result$mcmc[[2]][, 1:10], 
                                             FA.result$mcmc[[3]][, 1:10])) * 10)

#### estimate the model with varying number of topics from 60 to 80 (Fugire A.5) ####
for (i in 60:80) {
  # estimation results of the i-topic model
  load(paste0("STM_result_", i, ".Rdata"))
  ## estimate the Wordfish model for each topic
  topic.specific.positions <- matrix(NA, 10, i)
  rownames(topic.specific.positions) <- newspaper.names
  wordfish.result <- list()
  for (j in 1:i) {
    # document-feature matrix for topic j (taking a weighted sum of the frequency of term)
    weighted.dfm <- Study2.dfm * tcrossprod(STM.result$theta[, j], rep(1, ncol(Study2.dfm)))
    # compress the document-feature matrix at the newspaper level
    dfm.matrix <- matrix(NA, 10, ncol(Study2.dfm))
    rownames(dfm.matrix) <- newspaper.names
    colnames(dfm.matrix) <- colnames(Study2.dfm)
    for (k in 1:10) {
      dfm.matrix[k, ] <- round(colSums(weighted.dfm[Study2.dfm@docvars$newspaper == newspaper.names[k], ]))
    }
    # remove words not used by two or more newspapers
    dfm.matrix <- dfm.matrix[, -1 * which(colSums(dfm.matrix > 0) < 2)]
    dfm.matrix <- dfm.matrix[rowSums(dfm.matrix) > 0, ]
    # convert a matrix to a dfm object
    pseudo.text <- rep("", nrow(dfm.matrix))
    for (k in 1:nrow(dfm.matrix)) {
      for (l in 1:ncol(dfm.matrix)) {
        pseudo.text[k] <- paste(pseudo.text[k], 
                                paste(rep(colnames(dfm.matrix)[l], dfm.matrix[k, l]), collapse = " "))
      }
    }
    Wordfish.dfm <- dfm(pseudo.text, what = "fastestword")
    # estimate the Wordfish model
    set.seed(12345)
    wordfish.result[[j]] <- textmodel_wordfish(Wordfish.dfm, sparse = TRUE)
    # record topic-specific positions of the newspapers
    newspaper.ID <- match(rownames(dfm.matrix), newspaper.names)
    topic.specific.positions[newspaper.ID, j] <- wordfish.result[[j]]$theta
  }
  ## estimate the factor analysis model
  # function to create initial values
  initial.fun.preliminary <- function(seed1, seed2) {
    set.seed(seed1)
    b <- rnorm(i)
    x <- rnorm(10, 0, 0.01)
    p <- runif(i, 0.5, 1)
    list(beta = b, x = x, inv.phi = p, 
         .RNG.name = "base::Wichmann-Hill", .RNG.seed = seed2)
  }
  # create initial values
  initial.values.preliminary <- initial.fun.preliminary(12345, 67890)
  # estimate a factor analysis model using JAGS
  FA.result.preliminary <- run.jags(model = "factor_analysis_JAGS_preliminary.R", 
                                    monitor = c("x", "beta", "phi"), 
                                    data = list(y = topic.specific.positions, n = 10, m = i), 
                                    inits = initial.values.preliminary, 
                                    n.chains = 1, burnin = 5000, sample = 1000, modules = "glm", 
                                    summarise = FALSE, thin = 10, method = "parallel")
  ## fix the signs of factor loadings of some items so that Chunichi's ideal point becomes negative
  # whether Chunichi's ideal point was estimated to be negative
  Chunichi.score.sign.preliminary <- sign(mean(as.matrix(FA.result.preliminary$mcmc[, 3]))) == -1
  # posterior means of factor loadings
  post.loadings.preliminary <- colMeans(as.matrix(FA.result.preliminary$mcmc[, 11:(10 + i)]))
  # IDs of five topics whose factor loadings should be positive to make Chunichi's ideal point negative
  positive.items <- order(post.loadings.preliminary, decreasing = Chunichi.score.sign.preliminary)[1:5]
  # IDs of other topics
  unconstrained.items <- c(1:i)[-positive.items]
  # function to create initial values
  initial.fun <- function(seed1, seed2) {
    set.seed(seed1)
    b <- rnorm(i)
    b[positive.items] <- runif(5, 10, 15)
    x <- rnorm(10, 0, 0.01)
    p <- runif(i, 0.5, 1)
    list(beta = b, x = x, inv.phi = p,
         .RNG.name = "base::Wichmann-Hill", .RNG.seed = seed2)
  }
  # create initial values
  initial.values <- list()
  for (j in 1:3) {
    initial.values[[j]] <- initial.fun(100 * j, 1000 * j)
  }
  # estimate a factor analysis model using JAGS
  FA.result <- run.jags(model = "factor_analysis_JAGS.R", 
                        monitor = c("x", "beta", "phi"), 
                        data = list(y = topic.specific.positions, n = 10, m = i, 
                                    positive.items = positive.items, 
                                    unconstrained.items = unconstrained.items), 
                        inits = initial.values, 
                        n.chains = 3, burnin = 5000, sample = 5000, modules = "glm", 
                        summarise = FALSE, thin = 20, method = "parallel")
  # convergence check
  print(which(gelman.diag(FA.result$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1))
  # posterior draws of ideal points
  ideal.points.draws <- t(rbind(FA.result$mcmc[[1]][, 1:10], 
                                FA.result$mcmc[[2]][, 1:10], 
                                FA.result$mcmc[[3]][, 1:10])) * 10
  rownames(ideal.points.draws) <- newspaper.names
  ## summary statistics
  # point estimates (posterior means)
  ideal.points.mean <- rowMeans(ideal.points.draws)
  # 95% credible intervals (highest posterior density intervals)
  ideal.points.CI <- HPDinterval(as.mcmc(t(ideal.points.draws)))
  # ascending order of the point estimates
  ideal.points.order <- order(ideal.points.mean, decreasing = TRUE)
  Wordshoal.list <- list(ideal.points.mean = ideal.points.mean, 
                         ideal.points.CI = ideal.points.CI, 
                         ideal.points.order = ideal.points.order, 
                         wordfish.result = wordfish.result, 
                         topic.specific.positions = topic.specific.positions, 
                         FA.result = FA.result)
  save(Wordshoal.list,
       file = paste0("Wordshoal_result_", i, ".Rdata"))
  print(paste0("Wordshoal using ", i, "-topic model is finished at ", date(), "."))
}

#### compare the results (Figure A.5) ####
cairo_pdf("Figure_A5.pdf",
          width = 6, height = 7, pointsize = 8, family = "Helvetica")
layout(matrix(1:20, 4, 5, byrow = TRUE))
par(mar = c(2, 1, 2, 1))
for (i in c(60:64, 66:80)) {
  load(paste0("Wordshoal_result_", i, ".Rdata"))
  names(Wordshoal.list$ideal.points.mean) <- rank(Study2.ideal.points.mean)
  dotchart(Wordshoal.list$ideal.points.mean[Wordshoal.list$ideal.points.order], pch = 19, xlim = c(-1.1, 1.6), 
           main = i, xlab = "")
  segments(Wordshoal.list$ideal.points.CI[Wordshoal.list$ideal.points.order, 1], 1:10, 
           Wordshoal.list$ideal.points.CI[Wordshoal.list$ideal.points.order, 2], 1:10)
}
dev.off()

## intraclass correlation
ideal.points.mean.matrix <- matrix(NA, 10, 21)
for (i in 60:80) {
  load(paste0("Wordshoal_result_", i, ".Rdata"))
  ideal.points.mean.matrix[, 81 - i] <- Wordshoal.list$ideal.points.mean
}
round(ICC(ideal.points.mean.matrix)$results$ICC[6], 3)